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Through detailed comparisons between Embedded Atom Method (EAM) and first-principles cal- 
culations for Al, we find that EAM tends to fail when there are large electron density gradients 
present. We attribute the observed failures to the violation of the uniform density approximation 
(UDA) underlying EAM. To remedy the insufficiency of UDA, we propose a gradient-corrected 
EAM model which introduces gradient corrections to the embedding function in terms of exchange- 
correlation and kinetic energies. Based on the perturbation theory of "quasiatoms" and density 
functional theory, the new embedding function captures the essential physics missing in UDA, and 
paves the way for developing more transferable EAM potentials. With Voter-Chen EAM potential 
as an example, we show that the gradient corrections can significantly improve the transferability 
of the potential. 

PACS numbers: 31.15.xv, 61.50.Ah, 62.20.F- 



I. INTRODUCTION 

Atomistic simulations have become an increasingly 
powerful tool in materials research and a worthy partner 
of theory and experiment. Among the great many atom- 
istic models, the Embedded Atom Method (EAM) 1 ' 2 has 
emerged as one of the most successful and versatile ap- 
proaches, representing the mainstay of empirical atom- 
istic simulations. To date, EAM has been applied to 
a variety of material systems, such as liquids, metals 
and alloys, semiconductors, ceramics, polymers, nano- 
structures, and composite materials. Examples of prob- 
lems that EAM has studied include structure, energet- 
ics and dynamics of lattice defects, 3,4,5 elastic response 
and phonons, 6,7,8 fracture and plastic deformation, 9,10,11 
surface and surface growth, 12,13,14,15 thermodynamics 
properties, 1 ' 1 and phase transitions, 1 ' etc. The applica- 
tions of EAM simulations have been reviewed in Rcf. I 8. 
The success and popularity of EAM are a consequence of 
its sound theoretical foundation - the density functional 
theory (DFT) and its simple analytical expression. The 
former assures that the essential physics be captured by 
EAM and the latter endows EAM with excellent numer- 
ical efficiency, in par with pair potentials. 

Despite its great success, EAM suffers from a major 
deficiency - the lack of transferability. Most of EAM 
models are only reliable in regimes for which they were 
parameterized; beyond the regimes of parametrization, 
the reliability of EAM potentials quickly deteriorates. As 
a result, the predicability of EAM is often questionable in 
defect systems and in non-equilibrium conditions where 
relevant physical quantities are not known accurately a 
priori and hence not included in the parametrization of 
the potentials. As to all empirical models, the lack of 
transferability of EAM is an indication that some theo- 
retical approximations of EAM model are not generally 
valid. 

In this paper, we show that the lack of transferabil- 



ity of EAM is attributable to the uniform background 
density approximation of EAM embedding function. We 
find that EAM fails whenever there are large gradients 
of electron density in the system. We overcome the de- 
ficiency of the uniform density approximation (UDA) by 
proposing a density gradient corrected EAM model which 
incorporates the gradients of the valence electronic den- 
sity in the embedding function. Specifically, we intro- 
duce additional terms into the embedding function which 
correspond to the density gradient corrections to the 
exchange-correlation and kinetic energy contributions. 
Motivated by the Perdew-Burke-Ernzerhof (PBE) 19 Gen- 
eralized Gradient Approximation (GGA) of DFT and 
the perturbation theory of "quasiatoms" 20 , the present 
model applies to an inhomogencous background density 
and has the correct limiting behavior as the exact energy 
functions. As a consequence, it extends the applicability 
of EAM and paves the way for developing more transfer- 
able potentials. 



II. FAILURES OF THE UNIFORM DENSITY 
APPROXIMATION 

First, we demonstrate that the failure of EAM can be 
linked to the presence of large gradients of electron den- 
sity by comparing EAM with first-principles DFT cal- 
culations. We establish this fact in bulk Al for which 
EAM is supposed to work very well. Several excellent 
EAM potentials 21 ' 22 ' 23 exist and they are used for com- 
parisons. Both elastic properties and stacking fault en- 
ergy of Al are calculated. We compute the cohesive en- 
ergy per atom and the stress tensor as a function of the 
right Cauchy-Green deformation tensor CV, (i, j = 1,2,3) 
for a primitive unit cell of bulk Al. There are six inde- 
pendent C elements, with the diagonal and off-diagonal 
elements varying from -0.28 to 0.28 and -0.18 to 0.18, 
respectively, giving rise to different elastic deformations. 
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In order to provide precise comparisons, we utilize the 
sparse grid method 24 - a novel algorithm that allows us 
to represent a fine high-dimensional mesh very efficiently. 
Specifically, we sample 483,201 data points in sparse grids 
which correspond to a regular grid with 65 points in each 
of the six dimensions of Cy and 65 6 « 7.5 x 10 10 points in 
total. Moreover, by taking advantage of the underlying 
symmetry of the system, we further reduce the number 
of data points from 483,201 to 24.567 for which we carry 
out first-principles and EAM calculations. 

The first-principles DFT calculations are based on the 
plane-wave and Projector Augmented- Wave method 25 as 
implemented in the Vienna Ab initio Simulation Pack- 
age (VASP) 26 ' 27 . We use PBE-GGA with a high plane- 
wave cutoff energy of 360 eV to obtain reliable en- 
ergy and stress. The fc-points are sampled according to 
Monkhorst-Pack method with the fc-point spacing less 
than 0.0252 A -1 . A Gaussian smearing of 0.1 eV is em- 
ployed to speed up the convergence of the calculations. 
As for EAM calculations, we employ three widely used 
EAM potentials, developed by Ercolcssi and Adams 21 , 
Mishin ct al. 22 , and Voter and Chen 21 . The first two po- 
tentials were constructed by fitting to both experimental 
and first-principles data, while Voter-Chen potential was 
fitted only to experimental data. Although some EAM 
potentials 21 were constructed with different motivations 
from the original EAM model, they all use the UDA in 
effect. For the convenience of presentation, we define 
AV = and AQ = ^p*, where V (V ) and fi (fi„) 

are the volume and solid angle of the deformed (unde- 
formed) unit cell. The solid angle is defined relative to 
the basis vectors of the unit cell. AV and AO charac- 
terize the volumetric and non-volumetric deformation of 
the unit cell. 

In Figure 1, we present the difference in the cohe- 
sive energy (AE = Seam — -EVasp) and the stress tensor 
(||Aer|| = 1 1 cteam — cvaspID- Overall, we find excellent 
agreement between the first-principles and EAM results 
over a wide range of deformations - a remarkable feat 
of EAM. Furthermore, the errors are insensitive to the 
solid angle, Afl, which reflects the delocalized nature of 
the metallic bonds in Al, and thus justify the use of an 
angular - independent model in Al. On the other hand, 
we find that the EAM errors depend very sensitively on 
the change of volume, AV; in particular, the EAM val- 
ues deviate significantly from the first-principles results 
for large compressions. The errors in energy can reach as 
high as 1 eV/atom for 40% compression. This dramatic 
difference cannot be accounted for by the fitting errors of 
EAM because all three potentials show exactly the same 
behavior. Moreover, the shortest interatomic distance in 
the compressed unit cell is 2.1 A, which is still within the 
fitting range of the potentials. For example, the fitting 
range of bond length is from 2.0 A to 6.3 A in Mishin po- 
tential and 2.0 A to 5.6 A in Ercolessi-Adams potential. 
The results suggest that the errors come from the model 
itself. 

When the interatomic distance decreases, the gradient 
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FIG. 1: The difference for the elastic properties calculated 
by EAM and VASP. (a), (c), (e), and (g) show the cohesive 
energy differences in eV/atom; (b), (d), (f), and (h) show 
the stress differences in eV/A 3 . (a), (b), (c), and (d) are 
calculated by Ercolessi-Adams potential, (c) and (d) are the 
corresponding contour plots of (a) and (b), respectively, (e) 
and (f ) are calculated by Mishin potential and (g) and (h) are 
calculated by Voter-Chen potential, (e), (f), (g), and (h) are 
contour plots. 



of electronic density increases. For a large compression, 
the electron density gradient could become too large for 
the UDA of EAM to be valid. Indeed, we find that the 
density gradient increases considerably in compressions, 
with the maximum value of the gradients rising from 0.38 
A~ 4 for the perfect lattice to 0.54 A" 4 for 40% compres- 
sion. On the other hand, the expansion of the lattice 
reduces the density gradient, and thus does not violate 
the UDA. 

To further the argument, we perform additional cal- 
culations for the generalized stacking fault energy (7) 
surface, which along with elastic constants, determines 
the plastic behavior of materials. We have carried out 
293 energy calculations for the entire 7-surfacc with the 
sparse grid representation. The supercell consists of 9 
layers in the (111) direction for both EAM and VASP 
calculations. 

The 7-energy along [121] and [101] directions is shown 
in Figure 2a and 2b. Again, overall agreement between 
the three EAM potentials and VASP is good. However, 
in the neighborhood of the [101] unstable stacking fault 
and the run-on stacking fault (the last two entries in Ta- 
ble I), the magnitude of the energy error is significant. In 
particular, the largest error of EAM occurs at the run-on 
stacking fault in which the atoms in the two neighbor- 
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FIG. 2: 7-energy along (a) [121] and (b) [101] directions from 
VASP and EAM calculations. The horizontal axis is in the 
unit of the Burgers vector of Al (2.86 A). The contours of 
valence electronic density and the density gradient for the 
"run-on" configuration are shown in (c). The arrow represents 
the direction and magnitude of the density gradient. 



III. DENSITY GRADIENT CORRECTION 
MODEL 

Having established the importance of the density gra- 
dient, we propose a gradient corrected model which 
could potentially improve the transferability of EAM. 
The model is based on the pioneering work of Stott and 
Zaremba 2 " on "quasiatoms" . Stott et al. have shown 
that by using a perturbation expansion for an inhomo- 
geneous background density, the embedding energy of a 
"quasiatom" can be expressed rigorously as a function 
of the background density and its gradient. Based on 
the "quasiatoms" theory, we introduce three additional 
terms which account for the gradient corrections to the 
exchange, correlation and kinetic energy contributions to 
the embedding energy of EAM. In this context, the orig- 
inal embedding function of EAM can be regarded as the 
UDA to the embedding energy Specifically, the corrected 
embedding function becomes 



Fifaai) = F {p i ) + F c {p i )g{s i ) 

+F x {pi)h(s i )+F G (pi,s i ) > (1) 



ing (111) layers are right on top of each other, result- 
ing in large density gradients. The valence electronic 
density and its gradient in the "run-on" configuration is 
presented in Figure 2c. Noted that the maximum gradi- 
ent of valence electron density, |g| max of both the unsta- 
ble and the run-on stacking faults is comparable to the 
corresponding value of large compressions (<~ 0.5 A -4 ), 
suggesting that the failures of EAM can be indeed at- 
tributed to large density gradients, irrespective of the 
specific atomic configurations. Furthermore, the data in 
the brackets of the Table shows a general trend that the 
magnitude of EAM errors increases as the maximum den- 
sity gradient increases. Again the shortest interatomic 
distance in Table I is 2.33 A, which is within the fitting 
range of bond length for the EAM potentials. 



TABLE I: Fault vectors and energies for four stacking faults 
obtained from VASP and EAM calculations. All energies are 
in mj/m 2 . IqL denotes the maximum gradient of valence 
electronic density calculated by VASP, and is in A" 4 . PBE, 
EA, Mishin and Voter stand for the results calculated by 
VASP, Ercolessi- Adams, Mishin et al., and Voter-Chen EAM 
potential, respectively. The errors of EAM results are pre- 
sented in brackets. i? m m (A) is the nearest neighbor distance 
in the corresponding configurations. 



Vector PBE 



EA 



Mishin 



Voter 



Rn 



1/6[121] 111 121(10) 157(46) 87(-24) 0.398 2.86 

1/10[121] 184 132(-56) 190(6) 118(-66) 0.399 2.74 

1/4[101] 564 663(99) 603(39) 455(-109) 0.507 2.47 

1/3[121] 1208 1667(459) 1354(146) 1096(-112) 0.581 2.33 



where pi = pf" (Rij ) is the background density at atom 
i, and p^ is the density contribution from atom j. = 



Ri — Ra 



and Ri stands for the atomic coordinates. Al- 



though the background density pi is not the same as the 



total density p(r) where p (r) 



j 



(r~^l) 



they 



are closely related and the gradients of both densities are 
well defined. In particular, we can define a dimension- 



less background density gradient Si as Si tx 



/P- 



4/3 



where 



is the amplitude of background density 



gradient. In practice, Si can be approximated by its lo- 



cal average: Si ~ (s;) oc 



[p 



14/3 

1 jjti 



an, 



F (p) is the 



UDA embedding function and Fq is the gradient correc- 
tion to the kinetic energy. The leading term of Fq is of 
the von Weizsacker fornr x , and can be approximated as 
Fq (p, s) = A^o (p) k (s). Here Kq resembles the Thomas- 
Fermi kinetic energy 20. 



and k (s) = Xqj 



+ fe 2 lS ;i + fe22S 4 ' 

Ao, fcn, ki2, fc2i and AC22 are undetermined parameters. 

For exchange and correlation energy corrections, we 
adopt the functional form of PBE-GGA due to its sim- 
plicity. Fc and Fx in Equation 1 corresponds to the cor- 
relation and exchange energy of the local density approx- 
imation (LDA) of DFT and g (sj) and h (si) are the cor- 
responding gradient corrections. The explicit forms Fq 
and Fx can be found in standard references of LDA 31 ' 32 . 
We assume spin degeneracy here although the spin polar- 
ization can be considered easily and could be useful in the 
development of spin-dependent EAM potentials for mag- 
netic materials. In addition, we require that the modified 
embedding functions have the same limiting behavior as 
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^s 2 F x , 







-Fc + kqFx, s -> oc 



(2) 



and ft (s) = , which 



the exact functions: 
g(s)F c + h(s)F x 

We choose g (s) = - go+gi „ a+a<1 — w - 
satisfy the above conditions although other forms of g (s) 
and h (s) can also be used. 

Over all, there are six functions, Fc, Fx, Kq, p at (R), 
Fq (p), and ip (R), which are to be fitted. The first three 
are new terms and in conjunction with g(s), h(s) and 
k (s), they represent the gradient corrections to the em- 
bedding function. 

Since the UDA functions, Kq (p), Fc (p), and Fx (p), 
have the same functional forms as their DFT/LDA 



counterparts 



29,30,31,32 



we have 



k (p) cx p 5 / 3 . 



(3a) 



F C (p) =p(ci + c 2 r s ) In ( 1 + -t-^TT ) , r s oc 



M+ 1 
Fx(p)ccp 4 /\ 



1/3 ' 

(3b) 
(3c) 



Replacing the proportional sign 'oc' in all the above 
equations with an equal sign '=', one can fit the gradient 
corrected EAM (GCEAM) potential by introducing 13 
additional parameters. These additional parameters are 
ci, c 2 , 13, p, g , 9i, h , h x , A , k n , k 12 , k 21 , and k 22 . 
Among them, seven parameters (kn, k X2 , k 2 i, k 22 , go, gi 
and a) are introduced through the gradient corrections. 
kq and Aq can be absorbed into the embedding functions. 



IV. EXAMPLE: GRADIENT CORRECTED 
VOTER-CHEN POTENTIAL 

In this section, we apply the gradient corrections to the 
Voter-Chen (VC) potential, and the resultant potential is 
termed as GCEAM- VC. It is important to mention that 
the corrections are not constrained in any way by the 
specific form of the EAM potential. We choose the VC 
potential because its simplicity - it has only five param- 
eters, much fewer than other EAM potentials, such as 
Mishin and Ercolessi- Adams potentials. As a result, VC 
is not as accurate as the other EAM potentials. However, 
the simplicity of the VC potential renders more transpar- 
ent physics, and frees us from intensive parameter fitting 
- which is not the emphasis of the present paper. The goal 
of the article is to illustrate the importance of the den- 
sity gradient corrections in improving the transferability 
of EAM, rather than to generate the best possible EAM 
potential for Al. Had we started from an EAM potential 
with more parameters, we would have gotten even better 
test results for Al. Nevertheless, even with the VC po- 
tential, the gradient corrections can significantly improve 



the self-interstitial energies, stacking fault energies, etc. 
which involve high density gradient configurations. 

The GCEAM- VC potential takes the general form of 
EAM model. The cohesive energy of a system can be 
written as 



E = - ^2 p l3 (Rij) + ^2f, (pi, Si), 



(4) 



where embedding function Fi (pi,Sj) is expressed in Eq. 
(I)- 

The parameter fitting in GCEAM- VC follows the same 
procedure of the Voter-Chen potential , but with two 
modifications. The first modification is that the pairwisc 
interaction now is taken the form of: 



<p(R) 
<pi (R) 

V2 (R) 



Pi [R) + p 2 (R) . 



(5) 



D 



^ _ e -a M {R—R h 



D M , 



C 



R 



Here, epi (R) is a Morse potential used in the original 
Voter-Chen potential. ip 2 (R) is added to account for 
the repulsive interaction at short distance, (p (R — > 0) — > 
co when C > 0. However, if one prefers to use fewer 
parameters, p 2 (R) can be ignored without worsening the 
results (see the discussions of Fig. 4) . 

The second modification is that we do not fit the di- 
atomic molecular data. Instead, the force constants of 
bulk Al with different lattice constants are fitted because 
accurate force constants give rise to accurate phonon dis- 
persions, and hence, accurate thermal properties such as 
thermal capacity and conductivity. The force constants 
are fitted for several lattice constants, including 0.9ao, 
0.95ao, ao, 1.05ao, and l.lao, where eto is the equilibrium 
lattice constant of bulk Al. It is found the GCEAM- VC 
potential gives good description for the diatomic proper- 
ties without fitting them (See Table II). This is an exam- 
ple of improved transferability of the GCEAM model. 

The embedding function of both VC and GCEAM- VC 
potentials is determined by fitting the equation of states 
(EOS) to the universal EOS of Rose et al. 33 . Because 
the universal EOS does not agree exactly with the DFT 
(VASP) values (see Fig. 3), it is inevitable that both VC 
and GCEAM- VC potentials would deviate from DFT re- 
sults for large lattice expansions or large interatomic dis- 
tances. However, as shown later, the gradient corrections 
can improve significantly the description of high density 
gradient configurations involving lattice defects. 

The density function of the VC potential is given as 



p(r) =r 6 [e- &r + 2 9 e 



-2/3 2 r] 



(6) 



(3 2 needs to be fitted. We keep the same smoothness con- 
ditions for the pairwisc interaction, atomic density and 
EOS function in GCEAM- VC as in the VC potential 23 , 



5 




0.8 1 1.2 1.4 1.6 1.8 2 2.2 

a/aa 



FIG. 3: The cohesive energy per atom of fee Al as a function 
of the lattice parameter (the scaled equation of states) . E co h 
is the cohesive energy for equilibrium lattice constant ao- The 
Rose et al.'s universal EOS is represented by the blue line, and 
the VASP values are represented by red open circles. 

with the cutoff radius i? C ut of these functions to be fit- 
ted. Thus, there are seven parameters, Dm, olm, Rm, 
C, n, @2 and i? cu t to be fitted before applying the gra- 
dient corrections. Overall, there are twenty parameters 
in GCEAM-VC potential, including thirteen parameters 
associated with the gradient correction terms. The opti- 
mized values of all the parameters are presented in Ta- 
ble II. 

TABLE II: One set of optimized parameters for GCEAM-VC. 
Length and energy unit are in A and eV, respectively. 

ci -0.45618 gi -1.60375 k 12 -58.09505 R M 1.34426 

c 2 -1.03903 h 0.41370 k 2 i -10.15397 C 0.43750 

P 1.17295 fei 0.00327 fc 22 26.824 79 n 4.06038 

p 0.29559 Ao -0.52795 Dm 4.00963 fa 3.46742 

3o 1.63667 fen 157.54365 a M 2.05403 R cu t 5.56250 



The pair interaction tp(R) and the atomic density p(R) 
of GCEAM-VC potential are plotted in Fig. 4 in com- 
parison with the VC potential. It is found that the pair 
interaction of GCEAM-VC changes very little from that 
of VC. Although the atomic density function of GCEAM- 
VC potential appears to be rather different from that of 
VC, this turns out not to be the case. Using the fact that 
Eq. 4 is invariant under the transformation: 

p (R) -► tp (R) ,F(p,s)->F {pit, s/t) , 

we can define a scaled atomic density p (i?) = 
p (R) I max (p (R)). p (R) is plotted in the bottom panel 
of Fig. 4 and one finds little difference between the 
GCEAM-VC and VC atomic density functions. There- 
fore, we conclude that all improvements to the VC poten- 
tial come from the gradient corrections, i.e., the model 
itself. 

Some important properties predicted by VC and 
GCEAM-VC potentials are collected in Table III. From 
Table III, it is found that the GCEAM-VC improves 
the overall performance of the VC potential, especially 
for high density gradient configurations, such as sclf- 
intcrstitials. The table clearly demonstrates the success 
and improved transferability of GCEAM-VC. 
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FIG. 4: The pair interaction function tp(R) and atomic den- 
sity p(R) of VC (blue dash) and GCEAM-VC (red solid) po- 
tentials. 

Although GCEAM-VC gives a better result for the 
melting temperature than VC, the deviation from the ex- 
perimental value is still large. This is due to the fact that 
the melting process is associated with long-range inter- 
actions, whereas the density gradient corrections tend to 
be short-ranged. Therefore the gradient corrections are 
not expected to have significant effect on melting temper- 
ature. This is not an intrinsic problem of the GCEAM 
model because one could improve the melting tempera- 
ture by fitting more accurately the long-range tails of the 
UDA functions, e.g., <f>(R), p at (i?), and F (p). 

In Fig. 5, we compare the phonon dispersions between 
GCEAM-VC and VC, against the experimental data. It 
is found that GCEAM-VC predicts much better results 
than VC. 




FIG. 5: The phonon dispersion curves for Al. Red lines are 
calculated with GCEAM-VC potential, blue dash lines are 
calculated with VC potential, and open circles are experi- 
mental data taken from Ref. 40. 

Furthermore, with the help of the sparse-grid method, 
we calculate the entire 7-surfacc using first-principles 
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FIG. 6: The absolute error between the 7-surface from EAM potentials and first-principles calculations | Eq^ 1 — Eq^ 1 ' | . In 
(a), (b), (c) and (d), the EAM potentials are Ercolessi- Adams, Mishin, VC and GCEAM-VC, respectively. 



VASP method. From the 7-surface, one can derive the 
properties of all {lll}-type dislocations in Al 41 . The 
absolute errors between the 7-surface determined by var- 
ious EAM potentials and VASP calculations are shown 
in Fig. 6. The projection of the 7-surface along two 
special orientations is plotted in Figure 2a and 2b. It 
is found that the GCEAM-VC potential yields the most 
accurate result for overall 7-surface. Here we should em- 
phasize that GCEAM-VC does not fit any stacking fault 
configurations. In contrast, Ercolessi- Adams and Mishin 
et al. potentials both have included 7-cncrgies in their 
fitting database, and yet their results are not as good as 
GCEAM-VC. This is an important success of GCEAM 
potential in terms of transferability. Moreover, GCEAM- 
VC potential gives much more accurate stacking fault en- 
ergy near the "run-on" configuration ((-^jOJ, ^-^,1^, 

and f^p, 5^ in Fig. 6). These results confirm that in- 
deed the gradient corrections are crucial for describing 
high density gradient configurations, such as the "run- 
on" stacking faults. On the other hand, the errors of 
GCEAM-VC appear at configurations where interatomic 
distance is greater than that of a perfect lattice. These 
errors are not the intrinsic problem of the gradient cor- 
rection model, but rather due to the fitting strategy of 
the Voter-Chen potential. 

Finally, it is useful to mention that the force calcula- 
tion in GCEAM maintains the comparable numerical ef- 
ficiency with the standard EAM models. Thanks to the 



fact that the modified embedding functions can be fac- 
tored by a p-dependent term and an s-dependent term, 
the analytical expression of force remains simple - it has 
several additional terms that are of similar complexity 
of that of standard EAM. To compute these additional 
terms, the GCEAM needs to perform extra calculations 
of which the most time-consuming part is the second 
derivatives of the charge density with respect to distance. 
By using cubic spline interpolations, these calculations 
can be made rather efficient and as a result, the GCEAM 
force calculation takes less than twice of the CPU time 
of the standard EAM. The code package for calculat- 
ing the energy and force with GCEAM-VC potential is 
available via the World Wide Web 42 or via e-mail at wu- 
gaxp@gmail.com. 



V. DISCUSSION AND CONCLUSION 

Finally, it is instructive to relate the present correc- 
tions to other EAM models 2 ' 15 ' 44 . In the original EAM 
model, the electron correlations arising from the inho- 
mogeneous background density are largely ignored. The 
goal of the present model is to capture the missing corre- 
lations by taking into consideration of density gradients. 
Apart from the inhomogeneity of the density, the corre- 
lation effect also manifests itself in small molecules and 
clusters - a well-known fact in quantum chemistry that 
motivated the development of GGAs. By introducing 
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TABLE III: Properties of Al predicted by VC and and 
GCEAM-VC potentials in comparison with experimental 
and/or ab initio data. * Fitted properties. 



Experimental Voter-Chen 
or ab initio (Ref. 23) 



Lattice properties: 








fflo(A)* 


4.05 a 


4.05 


4.05 


£o(eV/atom)* 


-3.36* 


-3.36 


-3.36 


-B(GPa)* 


79 c 


79 


79.5 


cu(GPa)* 


114 c 


107 


113 


c 12 (GPa)* 


61.9 C 


65.2 


62.5 


c 44 (GPa)* 


31.6 C 


32.2 


32.9 


Diatomic Properties: 








£> e (eV) 


1.60 d 


1.54 


1.61 


Re{k) 


2.47 d 


2.45 


2.52 


Phonon frequencies: 








v L (X)(THz) 


9.69 e 


8.55 


9.62 


v T (A)(THz) 


5.80 c 


5.20 


5.36 


v L (L)(THz) 


9.69 e 


8.86 


10.1 


v T (i)(THz) 


4.19 c 


3.70 


3.82 


v L (if)(THz) 


7.59 c 


6.87 


7.69 


vry W(THz) 


5.64 c 


4.80 


5.00 


vt 2 (JC)(THz) 


8.65 e 


7.76 


8.69 


Vacancy: 








El{eW) 


0.68 / 


0.63 


0.65 


Self-interstitial: 








E{ (O h )(eV) 


2.73 


2.10 


2.41 


E[ (T d )(eV) 


3.08 


2.55 


2.87 


E{ ([111] dumbell)(eV) 


2.97 


2.48 


2.81 


E{ ([110] dumbell)(eV) 


2.76 


2.12 


2.38 


E[ ([100] dumbell)(eV) 


2.53 


2.02 


2.30 


Melting temperature: 








T m (K) 


933.6 


593.5±10 


672.5±: 



"Reference 34. 
b Reference 35. 
c Reference 36. 
d Reference 37. 
e Reference 38. 
■'Reference 39. 



a PBE GGA-like correction to the exchange-correlation 
part of the embedding energy, the present model could 
improve the description of the correlation effect. The 
modified EAM (MEAM) and its multistate variant strive 
to improve the transferability by making the background 
density p angular and reference-state dependent. How- 
ever, since they are based on UDA, the MEAM model 
does not treat the electron correlations adequately. As 
a result, it cannot deal with small clusters accurately 
as documented in the literature 4 '. The charge transfer 
EAM (CT-EAM) also recognizes the importance of the 
correlation effect. However it addresses the problem by 
introducing a reference-state (and its charge) dependent 
background density p. Since the present model consid- 
ers exchange-correlation energy explicitly, it can achieve 
the same goal of CT-EAM with a simpler function form. 
Moreover, one could incorporate MEAM and its variants 
into the present model by making the background density 
p in Eq. (1) as angular, reference-state and/or charge- 
dependent if so desired. 

In conclusion, we have performed detailed EAM and 
first-principles calculations of Al for elastic deformation 
and generalized stacking fault energy. We find that al- 
though EAM models reproduce well the first-principles 
results for most cases, they tend to fail when the elec- 
tron density gradients become substantial. We attribute 
the failures of EAM to the violation of UDA underlying 
the existing EAM models. To remedy the deficiency of 
UDA, we propose an improved EAM model which con- 
siders explicitly the gradient corrections to the embed- 
ding function in terms of the exchange-correlation energy 
and the kinetic energy. We show that the gradient cor- 
rected model can significantly improve the transferability 
of EAM, and represents a new direction for developing 
more transferable EAM potentials. 
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